PATENT APPLICATION 



TITLE OF THE INVENTION 
Long Elements Method for Simulation of Deformable Objects 



INVENTORS 
Remis Balaniuk 
Ivan F. Costa 
J. Kenneth Salisbury, Jr. 



BACKGROUND OF THE INVENTION 

1. Field of the Invention 

The present invention relates generally to deformable 
object modeling and, more particularly, to a novel long 
elements method (LEM) of real time simulation of deformable 
objects in a virtual environment. 

2. Description of the Related Art 

Deformable object modeling has been known in the art 
for decades and is widely utilized in fields such as 
engineering, computer-aided design (CAD) , and 
entertainment. Particularly in a virtual reality (VR) 
computing environment, the ability to model and manipulate 
deformable objects is essential to many applications. As 
such, graphic display of deformable objects has been 
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extensively studied in computer graphics and a number of 
methods have been proposed accordingly. These prior art 
methods ranging from non-physical methods, where individual 
or groups of control points or shape parameters are 
5 manually adjusted to shape editing and design, to methods 
based on continuum mechanics, which account for material 
properties and internal and external forces on object 
def ormation. 

As computational power increases, the ability to 
10 simulate physically based deformable models with enhanced 
multi-modal interactivity, i.e., using graphic and haptic 
I s * interfaces to manipulate and to change the topology of 
n deformable objects in real time, becomes even more 
'I" desirable. 

m 

;;jft 15 A major application area for this desirable real time 

physically based deformable object modeling with enhanced 
multi-modal interactivity is the simulation of 
biomaterials . Physical modeling of biomaterials has a 
broad range of applications ranging from understanding how 
20 soft and hard tissue respond under loading, to patient- 
specific planning of reconstructive procedures, to the 
training of surgical skills, and much more. Interaction 
with these models is necessary in order to impose 
conditions of interest, to examine results and alternative 
25 solutions, to learn surgical procedures by performing them 
in simulation . 

Given the intrinsic physical nature of these models it 
is desirable to have direct physical interactions with 
them. "Haptic" or direct physical interaction with 
30 simulated objects and/or subjects has become possible in 



'urn 



{'!§ 



S00-226/US 



2 



PATENT APPLICATION 



recent years because of advances in haptic (touch) 
interface technology, with devices now commercially 
available from companies such as the Immersion Corp. and 
SensAble Technologies. 

Not only does haptic technology enable convenient and 
expressive direct interactions with biomechanical 
simulations, it is also well suited to the intrinsically 
three-dimensional (3D) nature of biological models. 
Further, it links physical intuition to understanding of 
models, and, in training and rehearsal situations, it 
greatly enhances the learning process, much as in music 
training where one actually learns better and faster with 
the feeling and fingering of an instrument's keys. 

Modeling the biomechanics of muscles, tissues, and 
organs is intrinsically a computationally difficult 
undertaking - doing so at haptically real-time rates 
requires significant computational resources and 
algorithmic finesse . 

What is more, simulation methods must balance between 
two conflicting demands. Simulation for deformable objects 
can be classified by the degree of interaction they allow 
and their accuracy. The usefulness of a simulation method 
is defined by these two conflicting demands. Models 
focusing on interactiveness must have low latency and are 
based in some internal structure suitable for topological 
changes in real time. Models focusing on accuracy have the 
precision of their results limited only by some scale 
factor and the computational power. 

Unfortunately, these seemingly conflicting 
classifications frequently go together. Interactive 
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methods, such as the mass-spring method disclosed by G. 
Miller in "The motion dynamics of snake and worms 7 ', are 
mainly non-physical and often inaccurate, as pointed out by 
Zhuang et al. in "Haptic Interaction with Global 

5 Deformations", which is hereby incorporated herein by 

reference. Accurate physically based methods, such as the 
Finite Elements Method (FEM) , are typically simulated off- 
line and not in real time. What is more, modifying these 
methods to achieve real time performance usually 

10 compromises their accuracy and/or their interactiveness. 
Simulation methods can also be classified in 
accordance with additional aspects. For example, 

; ;5S | 

deformations can be dynamic or static, global or local. 

'l S5lf 

J* Models can be elastic or visco-elastic, linear or 

im 

lift 15 nonlinear, surface or volumetric. Collision detection and 

m 

12 handling can be implemented using different approaches. 
fi All these aspects define the way a simulation will behave, 

Q its accuracy and performance, and the kind of applications 

j JJ' for which it will be well suited. 

Q 20 The distinction between dynamic and static methods is 

of particular interest here. Dynamic methods simulate the 
evolving state of a physical system. The bodies have mass 
and energy distributed throughout. Differential equations 
and a finite state vector define each model. Numerical 

25 integration techniques approximate the system state 
(position and velocity) at discrete time steps. 

In static methods, such physical system is described 
by equilibrium equations or closed-form expressions, such 
as those taught by G.S. Chirikjian in "Closed-form 

30 primitives for generating volume preserving deformations", 
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which is hereby incorporated herein by reference. These 
static equations are solved to find a static solution at 
each time step. In static and quasi-static methods, the 
time and the system state are usually not considered. 
5 Thus, a state-based approach, i.e., a dynamic method, 

should be used for more comprehensive and accurate 
simulations of a deformable object/medium/subject, given 
significant internal dynamics of the deformable object, as 
well as the dynamics of the movements (translations and 
10 rotations) of the object in space. Note static and quasi- 
static methods are nonetheless useful in applications where 

U the objects deform but do not move, or move slowly, and the 

f!j deformable media is highly damped. 

I'* A number of schemes implementing these methods for 

I'fl 

lIH 15 deformable modeling have been developed, most notable ones 
£ are: mass-spring damper (MSD) systems, the boundary element 

» method (BEM) , the finite difference method (FDM) , and the 

i'i finite element method (FEM) , as described by Wu et al. in 

•; ssi? 

^ "Adaptive Nonlinear Finite Elements for Deformable body 

i B Qp 

□ 20 simulation using dynamic progressive meshes 7 ', which is 

hereby incorporated herein by reference. As is well known 
in the art and according to Wu et al., these prior art 
methods have either used approximate methods that are not 
physically accurate or linear methods that do not produce 

25 reasonable global behavior. 

For example, particle-based techniques are usually 
unstable and damping is extensively used to bring the 
system into a global equilibrium. In a particle-based 
simulation, damping and other constraints lead to a stiffer 

30 system, demanding shorter time steps to achieve stability. 
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The number and distribution of particles is also a source 
of problems. Larger meshes increase the stiffness of the 
system and become computationally very expensive. Smaller 
meshes increase inaccuracies and difficult to preserve 
volume. Uneven distribution of particles may easily 
generate unstable interaction forces and non-smooth 
graphical deformations. 

Finite element based methods generally produce 
globally accurate behavior but are usually too 
computationally expensive to be simulated in real-time. 
Pre-computation and condensation can be used to achieve 
real time rendering rates but changes in topology (re- 
meshing) demand an update of the stiffness matrix making 
the pre-computed data useless. 

To overcome these deficiencies, Wu et al. described an 
adaptive meshing scheme based on dynamic progressive meshes 
(DPM) using nonlinear finite elements. The DPM method has 
several limitations. For example, it relies on a hierarchy 
built offline. Thus, an initial detailed mesh is the limit 
that can be achieved in online refinement. Also, its tree 
structure has directionality and cannot split a node in any 
arbitrary direction locally. Consequently, the refined 
local region does not necessarily achieve optimal mesh 
quality. 

What is needed in the art is a new, computationally 
simple, fast, portable, and efficient simulation 
methodology capable of achieving real time performance and 
optimal mesh quality without compromising either accuracy 
or interactiveness, overcoming prior art drawbacks and 
limitations . 
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BRIEF SUMMARY OF THE INVENTION 
Accordingly, it is a primary object of the present 
invention to provide a novel long elements method (LEM) for 
5 real time physically based dynamic simulations of 

deformable media, the LEM is computationally simple, fast, 
portable, and efficient, capable of achieving real time 
performance with optimal mesh quality, interactiveness, and 
accuracy. 

10 Accordingly, it is also an object of the present 

invention to provide a new meshing strategy based on the 
M s inventive LEM, the LEM is well suited for real time 

n animation and virtual environment multi-modal interactions 
!|2 and particularly useful for soft tissue real time 

CR 15 simulation using graphic and haptic rendering. 
j|V It is another object of the present invention to 

!' provide a static solution based on the inventive LEM for 

Q elastic deformations of objects filled with uncompressible 

!«' fluid, which is a good approximation for biological 

O 20 tissues, wherein the static solution provides globally and 

U 

physically consistent deformations without pre-calculations 
or condensations. 

It is yet another object of the present invention to 
provide a novel combination of dynamic and static 
25 solutions, wherein the link between a static solution and a 
dynamic simulation is made using the duality between 
pressure (stress/strain) and force. 

It is therefore an object of the present invention to 
provide real time physically based modeling systems, 
30 methods, and apparatuses utilizing the inventive long 
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elements for dynamic elastic and plastic simulations of 
deformable media, objects, and subjects, the systems, 
methods, and apparatuses combining dynamic and static 
solutions, having a static state-less deformation problem 
5 solving means for calculating respective strains caused by 
internal and external stresses and for providing a static 
solution thereto, and having an computational means for 
integrating forces generated by the respective deformation 
to update the state respectively, 
10 It is a further object of the present invention to 

provide a static state-less deformation engine that is 
I'* based on a static solution for elastic deformations of 
O objects filled with uncompressible fluid, wherein the 

'.T volumes are discretised in a set of the inventive long 

CP 15 elements and an equilibrium equation is defined for each 
jV such element using dynamic variables including pressure, 

" density, volume, and stress, and wherein the set of static 

□ equations, plus Pascal's principle and volume conservation, 
!S are used in defining, solving, and finding object 

□ 20 deformations and forces, thereby obtaining globally and 

It siiif 

physically consistent deformations. 

Yet another object of the present invention to provide 
systems, methods, and apparatuses relating to defining an 
original and efficient three-dimensional (3D) meshing 

25 strategy based on the inventive long elements, the novel 
meshing strategy permitting approximation of state of and 
strain at any point on a volume with a reduced number of 
explicitly updated points, wherein number of elements per a 
meshed model is proportional to square length of a side of 

30 the model rather than its cube, thereby requiring 
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substantially less computational power and resources than 
cubes or tetrahedral based meshing strategies. 

It is therefore an object of the present invention to 
provide systems, methods, and apparatuses based on the 
inventive long elements, the systems, methods, and 
apparatuses integrating the simplicity and compliance of 
the static method to a full 3D dynamic simulation of 
deformable media, including localized and/or global 
deformations of elastic, plastic, and/or highly deformable 
material . 

Still further objects and advantages of the present 
invention will become apparent to one of ordinary skill in 
the art upon reading and understanding the following 
drawings and detailed description of the preferred 
embodiments. As it will be appreciated by one of ordinary 
skill in the art, the present invention may take various 
forms and may comprise various components and steps and 
arrangements thereof. Accordingly, the drawings are for 
purposes of illustrating a preferred embodiment (s) of the 
present invention and are not to be construed as limiting 
the present invention. 

BRIEF DESCRIPTION OF THE DRAWINGS 
Fig. 1 shows soft-tissue simulations according to an 
embodiment of the present invention. In these 
simulations, (a) bend; (b) palpation; and (c) 
twist, global and physically consistent 
deformation are obtained. 
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Fig. 2 shows a deformable medium, a cube, meshed in the 
inventive long elements according to the 
principles of the present invention. 

Fig. 3 shows modeling a deformable medium touched by a 
5 rigid probe according to an aspect of the present 

invention, wherein (a) shows the meshed model in 
wire mode; (b) shows three orthogonal planes 
crossing the deformable medium, each plane 
consisting of the inventive long elements; and 
io (c) shows one of the three planes that is 

affected by the rigid probe. 

* Fig. 4 shows another aspect of the present invention 

3 

i where (a) shows a deformable medium being pulled 

t by a probe or force and (b) shows a plane of long 

n 

n. 15 elements affected by the probe or force. 

Fig. 5 is a schematic drawing of the inventive mass-less 
long element according to an embodiment of the 
present invention . 
Fig. 6 illustrates the effect of surface tension in a 
20 simulation system according an aspect of the 

present invention . 
Fig. 7 illustrates yet another embodiment of the present 
invention where simulation of a deformable medium 
occurs in two different types of space 
25 simultaneously. 

Fig. 8 is a schematic drawing of the inventive spring- 
like dynamic long element according to another 
embodiment of the present invention. 
Fig. 9 shows exemplary LEM based deformation modeling. 

30 
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DETAILED DESCRIPTION OF THE INVENTION 
Known methods, procedures, systems, components, as 
well as computing environment, may be discussed or 
5 illustrated in the drawings, description and discussion 

herein without giving details, so as to avoid obscuring the 
principles of the invention. 

A survey of deformable modeling in computer graphics 
can be found in "A survey of deformable models in computer 
10 graphics'" by Gibson et al. Exemplary teachings on 

deformable modeling employing different prior art methods, 
i i * including the interactive methods, the "Geometric Nonlinear 
q Finite Element Method", and the "Boundary Element Method", 

>' w can be found in the following publications, which are 

!|! 15 hereby incorporated herein by reference: 

1. James et al., "Artdefo accurate real time deformable 
11 objects", Computer Graphics, v. 33, 65-72, 1999. 

H 2. Bro-Nielsen et al., "Real-time volumetric deformable 

HI- 

|;J- models for surgery simulation using finite elements and 
Q 20 condensation", Proceedings , Eurographics - Computer 
Graphics Forum, 57-66, 1996. 

3. Cavusoglu et al., "A lapraroscopic telesurgical 
workstation", IEEE Transactions on Robotics and Automation, 
15(4), 728-739, 1999. 
25 4. Vuskovic et al., "Realistic force feedback for virtual 
reality based diagnostic surgery simulators", Proceedings , 
IEEE Intl. Conf. On Robotics and Automation - ICRA, 1592- 
1598, 2000. 

5. Downes et al., "Virtual environments for training 
30 critical skills in laparoscopic surgery", Proc. Medicine 
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Meets Virtual Reality, 6, 316-322, 1998. 

6. d'Aulignac et al., "A haptic interface for a virtual 
exam of the human thigh", Proceedings , IEEE Intl. Conf. On 
Robotics and Automation - ICRA, 2452-2457, 2000. 
5 7. Cotin et al., "Real-time elastic deformations of soft 
tissues for surgery simulation", IEEE Transaction on 
Visualization on Computer Graphics, 5(1), 62-73, 1999. 

As previously discussed, deformations can be dynamic 
or static, global or local. Models can be elastic or 
10 visco-elastic, linear or nonlinear, surface or volumetric. 
Collision detection and handling can be implemented using 
i'* different approaches. All these aspects define the way a 
simulation will behave, its accuracy and performance, and 

n lad 

5* the kind of applications for which it will be well suited. 
ffi 15 The present invention is well suited for soft tissue 

i;h 

u , real time simulation, particularly for surgical simulation. 

11 Deterministic factors for such suitability include: 

□ unrestricted multi-modal interactiveness, including 

[ SJ interactive topological changes (cutting, suturing, 

l!3.20 removing material, etc.), physically based behavior, 
volumetric modeling (homogeneous and non-homogeneous 
materials) and scalability (high accuracy when needed) . 

According to an aspect of the present invention, a 
static solution based on the inventive long elements is 
25 provided for elastic deformations of objects filled with 
uncompressible fluid, which is a good approximation for 
biological tissues. The volumes are discretised in a set 
of static Long Elements (LE) , and an equilibrium equation 
is defined for each element using dynamic variables. These 
30 dynamic variables or "states" are used to collectively 
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refer to quantities that may vary significantly during 
simulation such as: pressure, volume, stress, strain, 
position, density, velocity, etc. The set of static 
equations, plus Pascal's principle and the volume 
conservation/preservation, are used to define a system that 
is solved to find the object deformations and forces. As a 
result, global and physically consistent deformations are 
obtained. Detailed mathematical solution, as well as a 
generic soft tissue VR simulator implemented in C++ in a 
Windows NT platform are disclosed in Applicants' recent 
publication, "LEM - An approach for physically based soft 
tissue simulation suitable for haptic interaction", 
Conference Paper, Fifth PHANTOM Users Group Workshop, 
October 2000, which is hereby expressly incorporated herein 
by reference . 

According to another aspect of the present invention, 
a combination of dynamic and static approaches is provided 
for dynamic simulation of deformable media using the 
inventive long elements. According to the teachings 
disclosed herein, a static state-less deformation engine in 
a state-based system solves the deformation problem for a 
given input of internal and external stresses calculating 
the respective strains while a dynamic simulation engine 
integrates the forces generated by the deformation to 
update the state of the system accordingly. 

To simulate a state based system according to another 
embodiment of the present invention, the novel static LEM 
is extended and linked to a dynamic method based on the 
same long elements meshing strategy. The present invention 
further provides a novel 3D integrator for providing a 3D 
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dynamic simulation and for linking the dynamic and static 
simulations . 

The present invention combines the simplicity and 
compliance of the static method with a 3D dynamic 
simulation of the deformable media. Highly deformable 
material, whether compressible or incompressible, can be 
simulated. Localized or global deformations can be 
simulated. The inventive long elements define an original 
and efficient 3D meshing strategy that permits the 
approximation of the state and the strain at any point in a 
volume from a reduced number of explicitly updated points 
while preserving/conserving the volume. The reduced, n 2 
complexity of a LEM based meshing strategy, where n is the 
length of a side of the meshing subject, greatly shortens 
corresponding computing time and thus enhances simulation 
speed and performance, allowing fast and realistic 
interactive applications. For example, the complexity of 
collision detection may be simplified utilizing the LEM. 
Similarly, a LEM based simulation can simplify making 
topology changes such as cutting and removing material, 
enabling real time feedback and making it well suited for 
virtual surgical procedures. Furthermore, simulation 
models based on the inventive LEM may be superimposed. 
Superposition of these models permits modeling 
inhomogeneous materials such as hard tumor in soft tissue. 
Potential LEM applications include industrial design, 
character and animated object simulation, real time 
simulation of objects ranging from furniture to human 
organs . 

Fig. 1 illustrates exemplary LEM based simulations, 
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showing the bend of an elastic rod (a) , a palpation of a 
soft object with a rigid haptic probe (b) , and a twist of 
the rod (c) . 

5 1. Long Elements (LE) 

The state of a deformable object at a given instant t 

is ultimately defined by the state of all its points at t. 

In order to simulate such an object in a computer at 

discrete time steps simplifications and approximations need 
10 to be done. In particle-based methods, the physics of a 

finite number of points of the object are simulated and a 

j.4 collective behavior is obtained with simplified links 

H 

between them (springs, dampers, etc.). 

M In finite element methods, limited volumetric forms 

,; 

j, :t X 

ijfl 15 inside the object are simulated (cubes, tetrahedral, etc.) 
; lM and a continuous behavior is obtained with boundary 

ii conditions at points shared by multiple elements. 

|ii£i 

.is* The LE method explores a different kind of 

'If?!? 

Ijj- approximation, based on the decomposition of the object 
Q 20 internal mechanics. This decomposition is based on three 
reference planes that cross the object and on the 
simulation of the relative positions of points inside the 
object with respect to these reference planes. 

25 1.1 Static Mesh 

Consider that at a given instant t in time three 
orthogonal planes P x , P y , and P z cross a deformable object 
or medium and suppose that a point 0 at the intersection of 
the 3 planes is inside the object such that the position of 

30 any point of the object can be uniquely defined by its 
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distances to the reference planes P x , P y , and P 2 . 

Define S as a set of points inside the object 
belonging to P x , P y , or P z at t. Suppose that the rest of 
the points inside the object are connected to the points in 
5 S by imaginary threads. Suppose that for each point P in S 
exist two threads connecting points between P and the 
surface of the object in opposite directions. A thread can 
be a straight-line segment starting in P and normal to its 
reference plane, for example, or it can have an arbitrary 
10 form, as discussed herein. Each point on the object 

belongs to exactly three threads, each one originating from 
j.4 a different orthogonal plane. 

n 

On a finite implementation, the cross sections defined 
M by the reference planes are divided in polygons and the 
i;fv 15 threads become long rods, i.e., Long Elements (LE) , 
!'^. crossing the object. Fig. 2 shows the meshing of a cube in 
ii LEs normal to each reference plane. The three orthogonal 



medium touched by a rigid probe according to an aspect of 
the present invention. In Fig. 3, (a) shows the meshed 
model in wire mode; (b) shows three orthogonal planes 
crossing the deformable medium, each plane consisting of 
25 the inventive LEs normal to the plane; and (c) shows one of 
the three planes that is affected by the rigid probe. As 
shown in Fig. 3, only LEs normal to the plane being 
affected by the rigid probe are deformed. Note although in 
Fig. 3 the reference planes are three mutually orthogonal 
30 planes, more or less than three planes may be used and they 

S00-226/US 16 




planes intersect at point 0. 

Figs. 3 and 4 show examplary meshing strategies based 
on the inventive LEs. Fig. 3 shows modeling a deformable 
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may or may not be orthogonal. Furthermore, the planes may 
be replaced in general by surfaces. The LEs may be used to 
represent objects of higher or lower dimension than three, 
i.e., planar objects, hypercubes, etc. In Fig. 4, (a) 
5 shows a deformable medium being pulled by a probe or force 
and (b) shows a plane of long elements affected by the 
probe or force. Again, only LEs normal to the plane being 
affected by the probe or force are deformed. 

10 1.2 Free Form Long Elements 

The simplified straight LEs on Fig. 2, normal to the 
|>4 reference plane, define a 1:1 relation between length of 

}Z the LE and the distance from its reference plane. 

■'i SI!? . 

M- ! A generic LE can be defined if for each point (or 

m 

Ijjff 15 section in a finite implementation) inside the LE the first 

I'fl 

order derivative dd=dL of the distance d with respect to 

j[ jsij? 

n the partial length L is known. In this case, the LE can 

!>«t have any arbitrary form. The approximation of the distance 

•I 

H 3 

^ from the plane at any point of a LE can be done by 

Q- 20 integrating Zdd/dLAL at discrete steps AL. 

In general, the LEM "meshes" or represents an object 
modeled in one space, e.g., 3D, by projecting the object's 
original representation (shape, state, properties, etc.) 
into a plurality of representations in lower dimensional 
25 space (s). Simulation of objects' behavior can be performed 
to advantage in these several lower dimensional spaces and 
the results are thereafter combined and projected back into 
the original higher dimensional space. As such, different 
meshing strategies can be conceived to fill/mesh/represent 
30 different objects with different forms of LEs. For 
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example, in a Cartesian implementation, the LEs are 
parallel to the Cartesian axis and cross the deformable 
medium. One mesh defines one face of the medium. Multiple 
meshes can then be combined to obtain a full 3D deformable 
5 medium. 

2. Deformation Engine 

To simulate a deformable elastic medium, each LE is 
modeled as a spring. A static method disclosed herein is 
10 used to estimate the length of each spring-like LE based on 
the pressures and tensions applied thereto. In a system 

|<* utilizing such static stateless deformation method, the LEs 

□ 

are springs only and have no mass. In one aspect of the 

?™ present invention, the simulated materials are assumed to 

CP 

Sip 15 be incompressible. For example, a deformable object maybe 

i : h 

! ■ simulated as a volume filled with uncompressible fluid. 

* This implicitly defines the material's Poisson' s ratio to 

if sK? 

h be 0.5. However, by means well known in the field of the 

est 

^ mechanics of materials, the LEM could be modified to 

»»d 

20 explicitly take account of other appropriate values of 
Poisson' s ratio to enable modeling and simulation of 
compressible materials. 

2.1 Pressure and Stress 

25 Consider the mass-less long elastic element 

illustrated in Fig. 5. The force F per unit of area A is 
defined as pressure: P = F/A. The force per area unit 
producing the deformation is also the stress. 

For small forces applied, the stress s in a material 

30 is usually linearly related to its deformation (its change 
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in length in the long elastic object) . Defining elasticity 
E as the variable relating stress s and the fractional 
change in length: AL/L, it is possible to write: 

s = EAL/L. (1) 
Since the stress s is related to the fractional change 
in length, the force can be related to the elongation AL in 
the well known form: 

F = kAL (2) 
where k = AE/L (3) 

Note that k is not constant, but it depends on the 
length L. Note also pressure may be made time varying to 
enable simulation of, for example, pulsing fluid flow 
internal to a body. 



2.2 Static Solution 

The static condition states that the forces, or 
pressures, in one sense have a correspondent of the same 
magnitude in the contrary sense on each point of the 
surface of the object, or P± n t = Pext- The external pressure 
Pext on the surface is affected by the atmospheric pressure 
and by the stress when an elongation exists, so 

Pext = Patm + EAL/L. (4) 

Note the surface tension also affects the external 
pressure, as described herein. 

Considering that the object is filled by fluid, the 
internal pressure (Pint) is formed by the pressure of the 
fluid (without gravity) and the effect of the gravity 
acceleration (g), such that 

Pint = P fluid + dgh (5) 
where h is the distance between the upper part of the 
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fluid and the point where the pressure is calculated. 

From the last three equations, a continuous equation 
can be obtained as 

EAL/L - AP = dgh (6) 
5 where AP = P fluid ~ Patm- 

Another external pressure to be considered comes from 
the contacts between the object and its environment. Two 
situations are considered: a soft contact and a rigid 
contact between the object and its environment. In a rigid 
10 contact, the elongation AL is defined by the penetration of 
the contact in order to make the surface follow a contact 
j<4 position (y) . Accordingly, the equation (6) can be 

,'!* rewritten for the elements where there is external rigid 
i" 4 contact as 

cn 

if! 15 AL = y (7) 

IT! 

In a soft contact,- it is not the position but the 

» force of the contact that is being considered. To obey the 

I* 

p action-reaction law, the force applied to the external 

! ! 1 1 

! !f contact and to the object must have the same magnitude. 

Q--20 This means that the external pressure P env applied by the 

contact must be equal to P env = Pext - Pint- At the points on 
the object surface, where external contacts exist, a term 
is added to the right side of the equation (6), such that 

EAL/L - AP = dgh + P env (8) 

25 

2.3 Simulating Long Elements 

A volume is filled (meshed) using LEs. An equilibrium 
equation is defined for each element based on the stated 
principles expressed in equations (6) -(8). Global 
30 constraints are added in order to obtain a global physical 
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behavior. 

A LE can be compared to a spring fixed in one 
extremity with the other extremity attached to a point in a 
deformable object surface. These LEs are relaxed when the 
object is not touched, i.e., not deformed. The spring 
constant K± depends on the LE's length L, as demonstrated in 
the equation (3) . 

The force due to each LE has a magnitude given by the 
equation (2) in the direction of the spring. A LE does not 
occupy real space and has no mass. In one embodiment, the 
real space inside the deformable medium is occupied by some 
incompressible fluid of density d. 

To make the connection between the elements, two 
border conditions are applied: 

1. Pascal's principle states that an external pressure 
applied to a fluid confined within a closed container is 
transmitted undiminished throughout the entire fluid. 
Mathematically: AP± = APj for any i and j. (9) 

2. The fluid is considered incompressible. This means that 
the volume conservation must be guaranteed when there is 
some external contact to the object. The volume 
dislocated by the contact will cause the dislocation of 
the entire surface. In other words, the variation of 
volume due to the elements touched by the contact have 
to be equal to the sum of the volume created by the 
dislocation of all untouched elements to ensure the 
volume conservation 

2jAAL = 0 (10) 
i=i 

where n is the total number of long elements. This 
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equation ensures the conservation of volume defined by 
the number of the LEs, implying the conservation of the 
object volume. 

The volumetric discretisation adopted by the static 
solution based on the inventive LEs has two main 
advantages: the number of LEs used to fill a deformable 
object is one dimension less than in a discretisation based 
on tetrahedric or cubic elements; the simulation for 
graphic and haptic feedback is thus greatly simplified and 
can be directly derived from the LEs. No intermediate 
geometric representation of the deformable object is 
needed. The use of static equations avoids problems 
related to numerical integration, ensuring stability for 
the simulation. Additionally, no pre-calculations or 
condensations are used, enabling real time topology 
changes. The static solution is particularly well suited 
to topological changes such as cutting because it 
intrinsically preserves volume. 

2.4 Surface Tension 

Deformations cause a change on the surface area of the 
deformable object, even if the object volume is kept 
constant. Fig. 6 illustrates the effect of such surface 
tension. This change creates forces on the surface 
generating a surface tension. One of the effects of these 
forces in a deformable object is to make the contours of 
the surface smoother. 

To reproduce these forces due to the change on the 
surface area, linear elastic connections are included 
between neighbor long elements coupling their changes in 
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length. A number of terms are added to the right side of 
the equation 4 corresponding to the neighborhood considered 
around element. These terms are of the form P = F/A = kx/A, 
where x is the difference between the deformations of a 
5 long element and its neighbor and k is a local spring 

constant. Thus, for a given LE i, the term relating its 
deformation to the deformation of its neighbor j is 

k ±j (ALi - ALj) / A± (11) 

10 2.5 Mathematical Solution 

Equations (8), (9) and (11) define the final equation 
for a long element i untouched, i.e., P en v = 0, or in soft 

AW. 

contact with the environment. Considering n neighbors (ji 

toj„), 

lift 15 (Ei/Li + nk s /A)ALi - k s {AL n + ... + AL jn ) /A - AP 

y. = digihi + Penvi + Porthoi (12) 

| ! where the superficial spring constant k s and A were chosen 

3 constant for all elements to make the notation easier. 

Note "superficial spring constant" is used here to 
3 20 represent the stiffness between adjacent endpoints of 

neighboring LEMs . Note also k s can be defined locally by 
the averaged elasticity of the specific elements. The term 
Porthoi represents an internal stress defined herein. 

The untouched LEs and the LEs in soft contact with the 
25 environment (12) plus the LEs in rigid contact with the 

environment (7) define a set of N equations, where N is the 
number of LEs used to fill the object. Adding the equation 
of volume conservation (10) , there are N + 1 equations and 
N + 1 unknowns: the pressure {AP) and the deformation of 
30 each element {AL± for i = 1 to N) . 
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These N +1 equations can be written as a problem of 
the type A.x = B where A is a sparse matrix. The problem 
and thus the corresponding simulation system can be solved 
using fast standard numerical methods. If rigid contacts 
5 are considered, the A matrix changes when the set of 

constrained LEs change, requiring a new matrix inversion 
A" 1 . However, if only soft contacts are considered, as is 
the case for typical medical simulation scenarios, the 
changes occur only on the B vector of pressures and gravity 
10 effects, and only matrix X vector multiplications are 
needed to obtain the deformation vector x. 

IIS 3. 3D Integrator 

ft «i 

According to an aspect of the present invention, 
m 15 simulation of a deformable medium occurs in two different 
types of spaces simultaneously: in one-dimensional (ID) LE 
space, and in three-dimensional (3D) Cartesian space. In 
the LE space, the LEs are always straight, parallel and 
homogeneous. As illustrated in Fig. 7, the 3D shape of an 

y 

20 object is derived from the configuration of the LEs . 

The transformation between the state of the LEs within 
the LE space to the state of the object in 3D is made 
considering one reference point 0 inside the object whose 
state in both spaces are coincident and estimating the 
25 position of all other points in 3D with respect to it. 

3.1 Dynamic Long Element 

The same massless LE described herein using a static 
formalism in equations (4) through (12) can be described as 
30 an energy-storing device. Elastic potential energy due to 
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either compression or stretching can be simulated to create 
forces and derive a state for each element. 

The LE is now modeled as a spring attached to a 
particle with a known mass. As shown in Fig. 8, A 

5 combination of two LEs attached to the same particle is 
used as the basic element to mesh the volumes. The LE is 
an one-dimension entity, defined by its length L. The 
particle attached to it is defined by a mass m and a state 
in an one-dimensional space as well. A deformed LE creates 

10 a force that is applied to the particle accelerating it. 
The mass of the particle corresponds to the mass of the 

l»* material inside the LEs connected to it. 

IJ 

q The particles belong to the set S of points as 

\* previously defined. Those are the only points where forces 

\ t J- 

i!p' 15 are known inside the object. For each particle, one scalar 

if j ■ 

y k value is given, representing a component of the force 

n vector at that point normal to the reference plane of which 

Q the particle belongs. An explicit integration of the state 

I'ij- 

! Sf of the particle is done at each time step in one dimension. 

SB!? 

p. 20 The result is two scalar values, displacement and velocity 

[r«Ss . 

of the particle with respect to the reference plane. 

According to the present invention, displacement means 

the length of an arbitrary trajectory and distance means 

the length of the shortest trajectory between two points 
25 (or between a point and a plane) . Furthermore, a cell is a 

unit of volume defined by 3 crossing LEs. Each point of 

the object belongs to only one cell. 

The state at any point of the object can be 

approximated using the state of its three crossing LEs. 
30 Each LE is used to estimate the distance from the point to 
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the LE reference plane. The estimation of this distance is 
made using the state of the LE 1 s particle plus the relative 
position of the point inside that LE, the partial lengths 
of the cells and the derivatives dd=dL of the distance d to 
the reference plane with respect to the length L. 

The 3D Integrator according to the present invention 
is capable of performing the following four tasks: 
integrating the ID state of the particles, estimating the 
distances from the particles to their reference planes (the 
particle 1 s position), estimating the numerical derivatives 
of the distance with respect to the length for each cell, 
and estimating internal stress at each cell. 

3.2 Generating Forces from Deformations 

On the LE space, only scalar values are being 
estimated: deformation, force, acceleration, velocity and 
position. The displacements of the particles are estimated 
at each time step. The dynamics on this space are very 
simplified. The LEs are modeled as springs attached to 
particles. For the dynamic LE configuration illustrated in 
Fig. 8, and using equations (2) and (3), the force F 
applied to a particle p attached to LEs L± and Lj is defined 
by 

F = k(AL± - ALj) . (13) 
This force is used to integrate the acceleration, velocity 
and position of the particle. The position corresponds to 
the displacement from the reference plane. 

3.3 Estimating the Position of a Particle 

The origin 0 belongs to all three reference planes and 
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is simulated using the same decomposition in three 
integrations. For each time step, the displacements of 0 
decomposed on the three directions and these displacement 
values correspond to the distances from 0 to the reference 
5 planes (this is true only for 0) . The reference planes are 
static and defined at the beginning of the simulation. 

For a point P of the object, if the distance | \0P\ I 
does not change between t and tl, and the velocity vectors 
P' + 0' then P described a spherical trajectory (a 
10 rotation) relative to 0 during this interval, that can be 
approximated by an arc connecting its initial and final 
j,4 relative positions. On the other hand, if the distance 
;,'S{ | \ OP\ | changes between t and tl, the trajectory can still 

M 5 be described as following a sphere having radius varying in 



20 difference between the displacements of 0 and P indicates a 

rotation of P about O. 

For each particle, two reference points are defined 

and they correspond to their respective closest point on 

each line defined by the intersection of its reference 
25 plane and each of the other two planes. We suppose that 

the distances between a particle and its reference points 

are always known. 

To estimate the relative position of a particle P with 

respect to 0 considering the spherical trajectory, compute 
30 the rotation of the particle about an axis crossing 0. A 



cn 
m 15 



time. 




To estimate the distance from a particle P to its 
reference plane, consider that the position of 0 indicates 
the translation of the whole object (inside the coordinate 
system defined by the reference planes) and that the 
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3D rotation of a point about an arbitrary axis is defined 
by 

r = M.q (14) 
where Mis a 3x3 matrix, q = [q x , q y , q z ] , the initial 
5 position, and r = [r x , r y , r z ] , the final position of the 
point . 

If a vector o = [o x , o y , o z ] represents the position of 
0 r then the final position of P considering translation and 
rotation is defined by P' = [p' x , p r yr p'z] = [r x + o xr r y + 
10 o yr r z + o z ] . M is defined as: 

M - uu T + (cos (a)) (I - uu T ) + sin(a)S (15) 
l«* where a is the rotation angle and u the rotation axis. 

The rotation axis u is defined by an unitary 

! : est 

directional vector u = v/||v|| = U' ,y' ,z') T . S is defined 
i:f! 15 as : 



FU- Considering i, j, e fx, y, zj and a particle P belonging 
q: to the plane P±, Rj is defined as the reference point at the 
intersection between plane P± and Pj and R k is defined as 
20 the reference point at the intersection between plane P± and 
P K . Two scalar values d jf d k are defined as dj = p - rj and 
dk = p - r k where p is the displacement of P, rj is the 
displacement of Rj and r k is the displacement of R k . The 
vector v is defined using v± = 0, Vj = d k and = -dj. 
25 The angle a of the rotation is defined by the 

difference p - o between the displacements of P and 0 and 
the distance d Pv between P and the rotation axis v. P is 
describing a circle around v, whose arc corresponds to p - 



0 -z' / 
z r 0 -x 1 
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O. 



From the equation of the arc of a circle: L 



where d Pv is the radius of the circle and L = p 



o, the 



5 



angle <j> is obtained: 

0 = (p - o) /d P „ 
Because P belongs to Pi then cfc = 0 . qj = 



Pi?j | | and 



(16) 



= | | PRjt | I . As a result, p'i gives the distance from the 
particle to the plane Pi. p' j and p' k are not relevant. 

The bend illustrated in Fig. 7 is a good example of 
the rotation of parts of the object. On the horizontal 
10 mesh space, from the center of the mesh to its extremities, 
a gradual increase of the relative displacements can be 
seen between vertical neighbors. In 3D, the result is a 
gradual inclination of the object caused by the increasing 
rotations of the LEs . 



3.4 Local Angles and Partial Derivatives 

The same rotational analysis is used at the cells 
level to determine internal deformations. The goal here is 
to determine cell by cell the deformation of a LE. A 



Q 20 length and a directional vector are defined for each cell 



in order to describe its shape. Corresponding partial 
derivatives are determined and internal stress on the cell 
is estimated. The change in length of each cell in a LE is 
defined by local pressure. Partial lengths are 

25 approximated using methods disclosed herein. 

The position of a cell c inside a LE corresponds to 
the sum of the partial lengths of all cells disposed 
between c and the particle. This position changes with the 
deformation of the LE. Neighbor LEs are those disposed 

30 side by side on the mesh. Neighbor cells are those 
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disposed at the same position in neighbor LEs when the 
object is at rest. 

Considering i, j, k 6 {x, y ,z} and a LE L± based on 
plane P±. Two neighbors L' ir L" \ are selected on the 
5 directions going from L± to its reference points. Suppose 
the directional vector at cell c ifn in Li is to be defined. 
In each neighbor LE, there is a corresponding neighbor cell 
c'i, n , c"i, n . The positions of the cells c± /nr c' ±, n , and c"i, n 
are p irDi p'±, n , and p"i, n , respectively. The displacements 
10 of the LEs L ir L' ±, L"± , i.e., displacement of the particles 
with respect to p± r are p if p r ±, and p" if respectively. Two 
!••* differences d lr d 2 are defined as: 

Q d lrB = {(p'i + P'i,n) ~ (Pi + Pi,n)) (17) 

!* d 2 ,n = {{P"i + P"i,n) ~ (Pi + Pi,n)) (18) 

lIH '15 The bend of Li at cell a, n can be described as a 

i'fi ■ 

| : V directional vector. The bending vector is approximated 
n using the same analogy of a rotation with respect to a 

p reference point. The differences d lfIit d 2 , n define a local 

! J; ' torsion of Li. The torsion is approximated as a rotation of 
Q 20 the cell about a point on the intersection between the 

three neighbor cells. Equation (14) is used to estimate 
the directional vector r n = r xr r y , r z . The input vector q 
is defined as: q± = 1; qj = 0; q k = 0 and the rotational 
axis v as v± = 0; Vj = d 2 , n , v k = -di, n . 
25 Using the same analogy of an arc of circle, angle a of 

the rotation is defined using the local "width" and 
"height" of the cell, the rotational axis and the 
differences d lrnr d 2 , n . At each cell n in L± there is a 
crossing with a LE Lj, n based on the plane P jf and with a LE 
30 L k/ri based in P k . There is a cell Cj, n , a at a position a 
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inside L jrn and a cell c k , n ,b at a position b inside L k , n 
corresponding to ci, n . Considering the "height" of cell c irn 
the length Ij, n , a of cell Cj, n ,a and the "width" of cell d, n 
the length h,n,b of cell Ck,n,b> 

a = ^n + d ln (19) 

The bending vector r is obtained using equation (14): 
<=cos (r 4±f-), 

j2 



r j = 



r k = — 



sin( + ) 
« n ( ) 



and 



aR 2 « + 4A* 

Distance from the reference plane to a cell is 
estimated as follows. The first order partial derivative 
ad/3/ at Ci, n is obtained directly from the bending vector: 
ad/a/ = n. The distance Pi, n from the reference plane P± to 
the outer face of cell c irn is estimated as: 



L j,a,o a 2,a l Ka^a 



a=l 



where / i/a is the length of cell c ifa . 

3.5 Internal Stress 

A cell is defined by its 3 local lengths, l x , l y , h, and 
its three bending vectors r x , r y/ r z that change during the 
simulation. Consider that the three bending vectors are 
supposed to define locally an orthonormal basis in space. 
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If three vectors r x/ r y , r z define the orthonormal basis, 
then the dot products: r x .r y = r x .r z = r y .r z = 0. 

Considering three crossing LEs: L K is based on plane 
P x , L y is based on plane P y , and L z is based on plane P z . 

5 The crossing happens on cells c x in L x , c y in L y and c z in 
L z . The lengths of these cells are: l x for c Xf l y for c yr and 
l z for c z . The bending vectors r x = {^,r/,r/} for c x , r y = {f^,r y z } 
for c y and r z = {r/,/f ,r z z } for c z are estimated. A vector of 
differences is defined as d = d ir dj,d k = { r x . r y/ r x . r z ,r y . r z } . 

10 A Jacobian matrix J is defined as: 



ddjl x ddjl y ddjl z 



ddj/l x ddj/l y dd s ll z 
dd k ll x dd k ll y dd k ll z 

j !P Using the inverse J" 1 of the Jacobian matrix J, a 

ft- vector of differences on the bending angles can be 

approximated: c = J" * d = { A / x , A / yr A / z } . 
|j 15 Using equation (1), the internal stress s± applied to 

ILj. cell Ci can be defined as: Si = G± * Mi/ h where Gi is the 

i s; t;. 

U shear modulus of cell c±. The total internal stress of a LE 
is thus the sum of the local stresses at each of its cells. 
The internal stress is represented in equation (12) as 

20 Portho • 

The bend illustrated in Fig. 7 further demonstrates a 
coordinated deformation of the meshes. On the horizontal 
mesh, there is a displacement of the external elements 
causing the rotation of the LEs, while on the vertical mesh 
25 a compression of the inner LEs defined smaller cells with 
respect to the outer neighbors bending the cells. The 
cells closer to the extremity have a bend more accentuated 
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in this case. 

3.6 Approximating Partial Lengths 

The meshing of the volume defines the initial 
dimensions of each cell. Physical modeling of the material 
defines the elasticity of each cell. At each time step 
during the simulation, a LE is submitted to external 
pressures from its contact with the environment, internal 
pressures of the liquid, internal pressures from the local 
stress at the cells level, superficial tension. 

The static simulation estimates one difference in 
length for the whole LE. This difference needs to be 
decomposed in differences at each cell. To approximate 
this decomposition, it is assumed that only the local 
stress varies from one cell to another and that all other 
sources of deformation act uniformly over the LE. 

The deformation on a LE composed by n cells caused by 
the local stress at the cells alone is defined as: 

d = % l l l /E l (21) 

i=l 

where Si is the local stress on cell i, h is the cell's 
length at rest, and E± is the elasticity of the cell. For a 
total deformation A/ of the LE, the length of cell i is 
estimated as: 

Mi = U + SihlEi + (A/ - d) In (22) 

4. System Implementations and Applications 

In a standard dual 700MHz PC, one iteration of the 
simulation loop takes about 0.05 seconds for a 600 elements 
mesh. Fig. 9 shows examples of deformations simulated 
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using the inventive LEM disclosed herein. 

In a system implemented with the teachings of the 
present invention, the global deformations are physically 
consistent and important phenomena such as the movement of 
all parts of the solid due to the preservation of volume 
are automatically produced. 

As in any dynamic simulation methodology based on 
explicit integration of differential equations, it is 
necessary to balance time steps and damping accordingly to 
the scales and spring factors used on the model in order to 
acquire stability. The present invention minimizes the 
stability problems and enables the simulation of stiffer 
models by providing a novel LEM for simulating a system 
defined by fewer particles than a standard 3D method. 

As will be understood by one of ordinary skill in the 
art, it is anticipated that the novel LEM may be 
implemented, updated and/or improved in various forms. For 
example, a formal analysis of the mapping between the LE 
space and 3D space can prove the correctness of the 
simulated geometry and physics. A study of the errors 
added by the first order approximations and linear elements 
is necessary to refine the LEM. A study of stability, ways 
to improve the LEM and its limits are also anticipated. 

The present invention can be implemented as part of a 
more general VR simulator. The static solution aspect of 
the LEM solves the deformation problem in an elegant and 
simple way, independently of other simulation aspects as 
movements, i.e., translation, rotation, dynamics 
integration, etc., topology changes and collision handling. 
To obtain a generic simulator, the LEM can be integrated to 
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methods suitable to handle these aspects. 

For example, to facilitate modeling and collision 
detection, the LEM can be coupled to implicit or parametric 
surface methods. The deformations of the model are 
5 simulated in a set of LE meshes defining deformable two- 
dimensional (2D) sheets and rendered as a 3D volume through 
an implicit of parametric formulation of the object shape. 
The transformation from 2D sheets to 3D surfaces can be 
done using texture mapping techniques. 
10 Of particular interest is the application of a 

surgical interface for performing virtual medical 
procedures. LEM based strategies, methodologies, and 
13 systems for meshing complex objects and algorithms for real 
|i time remeshing of objects to enable interactive topological 

IP.' is' changes (cutting, removing parts, suturing, etc.) are 
^ anticipated. 

f Furthermore, both elastic and plastic deformations can 

;3 be simulated using the principles of the present invention. 

5 To simulate an elastic medium, the base lengths of the 

itel • 

•3 20 element h are kept constant in equation (9) during the 
simulation. In a plastic deformation, the changes in 
length A/i change permanently the element: l±, n = h,n-i+M±, n -i 
where n indicates the time step of the simulation. 
Combinations of elastic and plastic deformations can be 
25 easily implemented as well. Note that for a plastic 
deformation, the system of static equations change, 
requiring a new matrix inversion. Using plastic 
deformation is also possible to mesh new objects. A 
"modeling clay'' can be implemented mapping points from an 
30 meshed object, such as a cube, to a new one, such as a 
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sphere. The differences in distance between the two sets 
of points are used to create external forces applied on the 
original object, molding it to fit the new shape. 

In a system implemented with the inventive long 
elements, the 3D internal state inside a deforming medium 
can be modeled at any point from a reduced number of 
explicitly updated elements. The reduction in one order of 
magnitude of the number of elements enables the simulation 
of more complex objects with less computational effort. 
The combination of a state less and a dynamic approach 
improves the compliance of the simulation. Large 
deformations that rapidly change the entire shape of the 
object can be simulated in a reduced number of time steps. 
Much of the simulation occurs on the LE ID space, where the 
complexity of all aspects of the simulation is strongly 
reduced. State-less deformations are estimated using only 
multiplications of matrices by vectors. The reduced number 
of elements further minimizes instability problems. 

Although the present invention and its advantages have 
been described in detail, it should be understood that 
various changes, substitutions, and alternations could be 
made herein without departing from the principles and the 
scopes of the invention. Accordingly, the scope of the 
present invention should be determined by the following 
claims and their legal eguivalents. 
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